##rm(list = ls()) 
install.packages("sandwich")
install.packages("lmtest")
install.packages("foreign") 
install.packages("data.table")
install.packages("car")
install.packages("readstata13")
install.packages("plotrix")
install.packages("descr")
install.packages("weights")
install.packages("miceadds") 
install.packages("mice")
install.packages("readxl")
install.packages("scales")
install.packages("ggplot2")
install.packages("gridExtra")
install.packages("grid")
install.packages("MASS")
install.packages("Hmisc")
install.packages("qte")
install.packages("cem")
install.packages("MatchIt")
install.packages("Zelig")
install.packages("RItools")
install.packages("dplyr")
install.packages("poliscidata")
install.packages("tidyr")
install.packages("sandwich")
install.packages("lmtest")




library(foreign) 
library(data.table)
library(car)
library(readstata13)
library(plotrix)
library(descr)
library(weights)
library(miceadds) 
library(mice)
library(readxl)
library(scales)
library(ggplot2)
library(gridExtra)
library(grid)
library(MASS)
library(Hmisc)
library(qte)
library(cem)
library(MatchIt)
library(Zelig)
library(RItools)
library(dplyr)
library(poliscidata)
library(lmtest)
library(tidyr)
library(sandwich)
library(lmtest)


# # Set the working directory to the location of the current script
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))  # Set working directory to the current script's directory
getwd()  # Check the current working directory

# load in file when it isnt in the environment already
load("HurricaneDataFinal2012to2014.RData")

HurricaneDataPlacebo = HurricaneData








#================================#
# TREATMENT COUNTY RELIEF POLICY #  
#================================#
# # select only the needed variables

dta_ReliefSurrounding <- HurricaneDataPlacebo
dta_ReliefSurrounding <- HurricaneDataPlacebo[, c("gender2", "Relief", "AgeCat2014", "Dem", "NPA", 
                                           "Third", "Black", "Hispanic", "Other", "EIPOnly2014", "EIPOnly2012",
                                           "VBMOnly2012", "VBMOnly2014", "EDOnly2012", "EDOnly2014", "AnyEarlyMethod2012", "AnyEarlyMethod2014",
                                           "Voted2012", "Voted2014", "CountyCode.x")]


# subset Relief Policy (Relief) (8 counties that received relief policies); control counties just the surrounding counties)

#============================================#
#++++++++++++++++++++++++++++++++++++++++++++#


#subsetting for 12 COUNTIES, only 8 Relief as treatment, and Limited to only 4 CONTROL counties surrounding

dta_ReliefSurrounding <- subset(dta_ReliefSurrounding, dta_ReliefSurrounding$CountyCode.x=="BAY" | dta_ReliefSurrounding$CountyCode.x=="FRA" | dta_ReliefSurrounding$CountyCode.x=="GUL" 
                                | dta_ReliefSurrounding$CountyCode.x=="JAC" | dta_ReliefSurrounding$CountyCode.x=="HOL" | dta_ReliefSurrounding$CountyCode.x=="WAK" 
                                | dta_ReliefSurrounding$CountyCode.x=="WAS" | dta_ReliefSurrounding$CountyCode.x=="CAL" | dta_ReliefSurrounding$CountyCode.x=="GAD"  
                                | dta_ReliefSurrounding$CountyCode.x=="LEO" | dta_ReliefSurrounding$CountyCode.x=="LIB" | dta_ReliefSurrounding$CountyCode.x=="WAL")

dta_ReliefSurrounding$Bay <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "BAY", 1, 0))
dta_ReliefSurrounding$Washington <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "WAS", 1, 0))
dta_ReliefSurrounding$Jackson <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "JAC", 1, 0))
dta_ReliefSurrounding$Calhoun <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "CAL", 1, 0))
dta_ReliefSurrounding$Gulf <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "GUL", 1, 0))
dta_ReliefSurrounding$Franklin <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "FRA", 1, 0))
dta_ReliefSurrounding$Gadsden <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "GAD", 1, 0))
dta_ReliefSurrounding$Liberty <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "LIB", 1, 0))
dta_ReliefSurrounding$Holmes <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "HOL", 1, 0))
dta_ReliefSurrounding$Wakulla <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "WAK", 1, 0))
dta_ReliefSurrounding$Leon <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "LEO", 1, 0))
dta_ReliefSurrounding$Walton <- with(dta_ReliefSurrounding, ifelse(CountyCode.x == "WAL", 1, 0))

table(dta_ReliefSurrounding$Relief)



#==========================================================#
# Difference in Difference Models for the Relief Policy counties and 4 surrounding counties #
#==========================================================#

# Voted2012 vs. Voted 2014 (any method, early or election day)

# Reshape dataset from wide to long
dta_ReliefSurrounding_long <- reshape(dta_ReliefSurrounding,
                                      varying = c("Voted2012", "Voted2014"),
                                      v.names = "voted",
                                      timevar = "year",
                                      times = c("Voted2012", "Voted2014"),
                                      new.row.names = 1:946286,    #multiply the number of obs by 2
                                      direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurrounding_long$year2014 <- ifelse(dta_ReliefSurrounding_long$year=='Voted2014', 2014, 2012)
dta_ReliefSurrounding_long$d14 <- ifelse(dta_ReliefSurrounding_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurrounding_long <- data.table(dta_ReliefSurrounding_long)


# DiD variable: treatment*year
dta_ReliefSurrounding_long$did <- dta_ReliefSurrounding_long$Relief * dta_ReliefSurrounding_long$d14

#table(dta_ReliefSurrounding_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurrounding_long)
# 
# summary(m_ReliefSurrounding)

m_ReliefSurrounding <- coeftest(m_ReliefSurrounding, vcov=vcovHC(m_ReliefSurrounding, type="HC0", cluster="CountyCode.x")) 

m_ReliefSurrounding

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurrounding_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.9674 -0.4628  0.1628  0.3377  0.8609 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6491697  0.0015084  430.364  < 2e-16 ***
#   d14         -0.1122619  0.0012302  -91.257  < 2e-16 ***
#   Relief      -0.0567112  0.0012988  -43.664  < 2e-16 ***
#   did         -0.0722536  0.0018273  -39.541  < 2e-16 ***
#   AgeCat20142  0.0933976  0.0013996   66.731  < 2e-16 ***
#   AgeCat20143  0.2560783  0.0013259  193.129  < 2e-16 ***
#   AgeCat20144  0.3098964  0.0014351  215.939  < 2e-16 ***
#   gender2     -0.0235947  0.0009196  -25.656  < 2e-16 ***
#   Dem         -0.0473426  0.0011234  -42.144  < 2e-16 ***
#   NPA         -0.1947494  0.0013538 -143.850  < 2e-16 ***
#   Third       -0.0744080  0.0064010  -11.624  < 2e-16 ***
#   Black        0.0083583  0.0013076    6.392 1.64e-10 ***
#   Hispanic    -0.0207154  0.0028978   -7.149 8.77e-13 ***
#   Other       -0.0505346  0.0021415  -23.598  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.441 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.127,	Adjusted R-squared:  0.127 
# F-statistic: 1.052e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# i.e. those living in counties hit by michael who voted in 2012 (by any method) were less likely to vote (by any method) in 2014 than were those living in counties not hit by michael who voted in 2012 (by any method)


### VOTING METHODS ###


# Voted 2012 vs. Any Early Method 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted12VsAnyEarly_long <- reshape(dta_ReliefSurrounding,
                                                       varying = c("Voted2012", "AnyEarlyMethod2014"),
                                                       v.names = "voted",
                                                       timevar = "year",
                                                       times = c("Voted2012", "AnyEarlyMethod2014"),
                                                       new.row.names = 1:946286,    #multiply the number of obs by 2
                                                       direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVoted12VsAnyEarly_long$year2014 <- ifelse(dta_ReliefSurroundingVoted12VsAnyEarly_long$year=='AnyEarlyMethod2014', 2014, 2012)
dta_ReliefSurroundingVoted12VsAnyEarly_long$d14 <- ifelse(dta_ReliefSurroundingVoted12VsAnyEarly_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted12VsAnyEarly_long <- data.table(dta_ReliefSurroundingVoted12VsAnyEarly_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted12VsAnyEarly_long$did <- dta_ReliefSurroundingVoted12VsAnyEarly_long$Relief * dta_ReliefSurroundingVoted12VsAnyEarly_long$d14

table(dta_ReliefSurroundingVoted12VsAnyEarly_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted12VsAnyEarlySurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVoted12VsAnyEarly_long)
# 
# summary(m_ReliefVoted12VsAnyEarlySurrounding)

m_ReliefVoted12VsAnyEarlySurrounding <- coeftest(m_ReliefVoted12VsAnyEarlySurrounding, vcov=vcovHC(m_ReliefVoted12VsAnyEarlySurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVoted12VsAnyEarly_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.9734 -0.4198  0.1096  0.3854  0.9476 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6593099  0.0015131  435.727  < 2e-16 ***
#   d14         -0.3730352  0.0012340 -302.294  < 2e-16 ***
#   Relief      -0.0543818  0.0013029  -41.740  < 2e-16 ***
#   did          0.0429501  0.0018330   23.431  < 2e-16 ***
#   AgeCat20142  0.0666737  0.0014040   47.489  < 2e-16 ***
#   AgeCat20143  0.2202640  0.0013301  165.602  < 2e-16 ***
#   AgeCat20144  0.3087043  0.0014396  214.439  < 2e-16 ***
#   gender2     -0.0232394  0.0009225  -25.192  < 2e-16 ***
#   Dem         -0.0446825  0.0011269  -39.652  < 2e-16 ***
#   NPA         -0.1637920  0.0013581 -120.607  < 2e-16 ***
#   Third       -0.0649231  0.0064210  -10.111  < 2e-16 ***
#   Black        0.0053499  0.0013117    4.079 4.53e-05 ***
#   Hispanic    -0.0081587  0.0029068   -2.807    0.005 ** 
#   Other       -0.0354116  0.0021482  -16.485  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4424 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.2055,	Adjusted R-squared:  0.2055 
# F-statistic: 1.871e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# i.e. those living in counties hit by michael who voted (by any method) in 2012 were MORE likely to vote early (by any method) in 2014 than were those living in counties not hit by michael who voted (by any method) in 2012 


# Voted 2012 vs. Election Day 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted12VsED_long <- reshape(dta_ReliefSurrounding,
                                                 varying = c("Voted2012", "EDOnly2014"),
                                                 v.names = "voted",
                                                 timevar = "year",
                                                 times = c("Voted2012", "EDOnly2014"),
                                                 new.row.names = 1:946286,    #multiply the number of obs by 2
                                                 direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVoted12VsED_long$year2014 <- ifelse(dta_ReliefSurroundingVoted12VsED_long$year=='EDOnly2014', 2014, 2012)
dta_ReliefSurroundingVoted12VsED_long$d14 <- ifelse(dta_ReliefSurroundingVoted12VsED_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted12VsED_long <- data.table(dta_ReliefSurroundingVoted12VsED_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted12VsED_long$did <- dta_ReliefSurroundingVoted12VsED_long$Relief * dta_ReliefSurroundingVoted12VsED_long$d14

table(dta_ReliefSurroundingVoted12VsED_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted12VsEDSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVoted12VsED_long)

# summary(m_ReliefVoted12VsEDSurrounding)

m_ReliefVoted12VsEDSurrounding <- coeftest(m_ReliefVoted12VsEDSurrounding, vcov=vcovHC(m_ReliefVoted12VsEDSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVoted12VsED_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.85065 -0.24089 -0.05102  0.24841  1.09366 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.7073516  0.0014168  499.262  < 2e-16 ***
#   d14         -0.4871646  0.0011555 -421.623  < 2e-16 ***
#   Relief      -0.0405730  0.0012199  -33.259  < 2e-16 ***
#   did         -0.0922082  0.0017163  -53.724  < 2e-16 ***
#   AgeCat20142  0.0640476  0.0013146   48.720  < 2e-16 ***
#   AgeCat20143  0.1433016  0.0012454  115.064  < 2e-16 ***
#   AgeCat20144  0.1313528  0.0013479   97.447  < 2e-16 ***
#   gender2     -0.0174272  0.0008638  -20.176  < 2e-16 ***
#   Dem         -0.0269462  0.0010551  -25.539  < 2e-16 ***
#   NPA         -0.1225683  0.0012716  -96.389  < 2e-16 ***
#   Third       -0.0470251  0.0060122   -7.822 5.22e-15 ***
#   Black       -0.0094365  0.0012282   -7.683 1.55e-14 ***
#   Hispanic    -0.0163964  0.0027218   -6.024 1.70e-09 ***
#   Other       -0.0410678  0.0020114  -20.418  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4142 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.3117,	Adjusted R-squared:  0.3117 
# F-statistic: 3.275e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted (by any method) in 2012 were less likely to vote on Election Day in 2014 than were those living in counties not hit by michael who voted (by any method) in 2012

# Voted 2012 vs. VBM Only 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted12VsVBM_long <- reshape(dta_ReliefSurrounding,
                                                  varying = c("Voted2012", "VBMOnly2014"),
                                                  v.names = "voted",
                                                  timevar = "year",
                                                  times = c("Voted2012", "VBMOnly2014"),
                                                  new.row.names = 1:946286,    #multiply the number of obs by 2
                                                  direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVoted12VsVBM_long$year2014 <- ifelse(dta_ReliefSurroundingVoted12VsVBM_long$year=='VBMOnly2014', 2014, 2012)
dta_ReliefSurroundingVoted12VsVBM_long$d14 <- ifelse(dta_ReliefSurroundingVoted12VsVBM_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted12VsVBM_long <- data.table(dta_ReliefSurroundingVoted12VsVBM_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted12VsVBM_long$did <- dta_ReliefSurroundingVoted12VsVBM_long$Relief * dta_ReliefSurroundingVoted12VsVBM_long$d14

table(dta_ReliefSurroundingVoted12VsVBM_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted12VsVBMSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVoted12VsVBM_long)

# summary(m_ReliefVoted12VsVBMSurrounding)

m_ReliefVoted12VsVBMSurrounding <- coeftest(m_ReliefVoted12VsVBMSurrounding, vcov=vcovHC(m_ReliefVoted12VsVBMSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVoted12VsVBM_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.91535 -0.19264 -0.01512  0.23714  1.12843 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6931649  0.0012839  539.884  < 2e-16 ***
#   d14         -0.6142461  0.0010471 -586.626  < 2e-16 ***
#   Relief      -0.0446918  0.0011055  -40.426  < 2e-16 ***
#   did         -0.0014815  0.0015554   -0.953   0.3408    
# AgeCat20142  0.0452023  0.0011913   37.944  < 2e-16 ***
#   AgeCat20143  0.1383371  0.0011286  122.574  < 2e-16 ***
#   AgeCat20144  0.2221850  0.0012215  181.893  < 2e-16 ***
#   gender2     -0.0246202  0.0007828  -31.453  < 2e-16 ***
#   Dem         -0.0239508  0.0009562  -25.049  < 2e-16 ***
#   NPA         -0.1094763  0.0011523  -95.003  < 2e-16 ***
#   Third       -0.0387170  0.0054483   -7.106 1.19e-12 ***
#   Black       -0.0142527  0.0011130  -12.805  < 2e-16 ***
#   Hispanic    -0.0058276  0.0024665   -2.363   0.0181 *  
#   Other       -0.0270825  0.0018227  -14.858  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3754 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.4252,	Adjusted R-squared:  0.4252 
# F-statistic: 5.349e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted (by any method) in 2012 were MORE likely to vote by mail in 2014 than were those living in counties not hit by michael who voted (by any method) in 2012

# Voted 2012 vs. EIP Only 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted2012VsEIP_long <- reshape(dta_ReliefSurrounding,
                                                    varying = c("Voted2012", "EIPOnly2014"),
                                                    v.names = "voted",
                                                    timevar = "year",
                                                    times = c("Voted2012", "EIPOnly2014"),
                                                    new.row.names = 1:946286,    #multiply the number of obs by 2
                                                    direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVoted2012VsEIP_long$year2014 <- ifelse(dta_ReliefSurroundingVoted2012VsEIP_long$year=='EIPOnly2014', 2014, 2012)
dta_ReliefSurroundingVoted2012VsEIP_long$d14 <- ifelse(dta_ReliefSurroundingVoted2012VsEIP_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted2012VsEIP_long <- data.table(dta_ReliefSurroundingVoted2012VsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted2012VsEIP_long$did <- dta_ReliefSurroundingVoted2012VsEIP_long$Relief * dta_ReliefSurroundingVoted2012VsEIP_long$d14

table(dta_ReliefSurroundingVoted2012VsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted2012VsEIPSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVoted2012VsEIP_long)

# summary(m_ReliefVoted2012VsEIPSurrounding)

m_ReliefVoted2012VsEIPSurrounding <- coeftest(m_ReliefVoted2012VsEIPSurrounding, vcov=vcovHC(m_ReliefVoted2012VsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVoted2012VsEIP_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.90747 -0.34080  0.01898  0.30259  1.01898 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6836367  0.0014567  469.306  < 2e-16 ***
#   d14         -0.5067270  0.0011880 -426.539  < 2e-16 ***
#   Relief      -0.0479336  0.0012543  -38.216  < 2e-16 ***
#   did          0.0674272  0.0017647   38.210  < 2e-16 ***
#   AgeCat20142  0.0587951  0.0013516   43.499  < 2e-16 ***
#   AgeCat20143  0.1894142  0.0012805  147.924  < 2e-16 ***
#   AgeCat20144  0.2166800  0.0013859  156.346  < 2e-16 ***
#   gender2     -0.0156912  0.0008881  -17.668  < 2e-16 ***
#   Dem         -0.0450177  0.0010848  -41.497  < 2e-16 ***
#   NPA         -0.1459265  0.0013074 -111.614  < 2e-16 ***
#   Third       -0.0637463  0.0061816  -10.312  < 2e-16 ***
#   Black        0.0071576  0.0012628    5.668 1.45e-08 ***
#   Hispanic    -0.0061709  0.0027984   -2.205   0.0274 *  
#   Other       -0.0342739  0.0020680  -16.573  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4259 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.2745,	Adjusted R-squared:  0.2745 
# F-statistic: 2.736e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted (by any method) in 2012 were MORE likely to early-in-person in 2014 than were those living in counties not hit by michael who voted (by any method) in 2012


# Election Day 2012 vs. EIP Only 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingEDOnly2012VsEIP_long <- reshape(dta_ReliefSurrounding,
                                                     varying = c("EDOnly2012", "EIPOnly2014"),
                                                     v.names = "voted",
                                                     timevar = "year",
                                                     times = c("EDOnly2012", "EIPOnly2014"),
                                                     new.row.names = 1:946286,    #multiply the number of obs by 2
                                                     direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingEDOnly2012VsEIP_long$year2014 <- ifelse(dta_ReliefSurroundingEDOnly2012VsEIP_long$year=='EIPOnly2014', 2014, 2012)
dta_ReliefSurroundingEDOnly2012VsEIP_long$d14 <- ifelse(dta_ReliefSurroundingEDOnly2012VsEIP_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingEDOnly2012VsEIP_long <- data.table(dta_ReliefSurroundingEDOnly2012VsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingEDOnly2012VsEIP_long$did <- dta_ReliefSurroundingEDOnly2012VsEIP_long$Relief * dta_ReliefSurroundingEDOnly2012VsEIP_long$d14

table(dta_ReliefSurroundingEDOnly2012VsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefEDOnly2012VsEIPSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingEDOnly2012VsEIP_long)

# summary(m_ReliefEDOnly2012VsEIPSurrounding)

m_ReliefEDOnly2012VsEIPSurrounding <- coeftest(m_ReliefEDOnly2012VsEIPSurrounding, vcov=vcovHC(m_ReliefEDOnly2012VsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingEDOnly2012VsEIP_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.3568 -0.2845 -0.2361  0.6537  0.8900 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.2552468  0.0014941 170.836  < 2e-16 ***
#   d14         -0.0368682  0.0012185 -30.257  < 2e-16 ***
#   Relief      -0.0441624  0.0012865 -34.328  < 2e-16 ***
#   did          0.0786164  0.0018100  43.435  < 2e-16 ***
#   AgeCat20142  0.0325878  0.0013863  23.506  < 2e-16 ***
#   AgeCat20143  0.0918677  0.0013134  69.949  < 2e-16 ***
#   AgeCat20144  0.0556894  0.0014215  39.177  < 2e-16 ***
#   gender2      0.0015942  0.0009109   1.750 0.080092 .  
# Dem         -0.0264095  0.0011127 -23.735  < 2e-16 ***
#   NPA         -0.0759330  0.0013410 -56.625  < 2e-16 ***
#   Third       -0.0341009  0.0063403  -5.378 7.51e-08 ***
#   Black        0.0080992  0.0012952   6.253 4.02e-10 ***
#   Hispanic    -0.0095858  0.0028703  -3.340 0.000839 ***
#   Other       -0.0251378  0.0021211 -11.851  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4368 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.01312,	Adjusted R-squared:  0.01311 
# F-statistic: 961.3 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted on election day in 2012 were MORE likely to vote early-in-person in 2014 than were those living in counties not hit by michael who voted on election day in 2012


# VBM Only 2012 vs. Voted 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVBMVsVoted_long <- reshape(dta_ReliefSurrounding,
                                                varying = c("VBMOnly2012", "Voted2014"),
                                                v.names = "voted",
                                                timevar = "year",
                                                times = c("VBMOnly2012", "Voted2014"),
                                                new.row.names = 1:946286,    #multiply the number of obs by 2
                                                direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVBMVsVoted_long$year2014 <- ifelse(dta_ReliefSurroundingVBMVsVoted_long$year=='Voted2014', 2014, 2012)
dta_ReliefSurroundingVBMVsVoted_long$d14 <- ifelse(dta_ReliefSurroundingVBMVsVoted_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVBMVsVoted_long <- data.table(dta_ReliefSurroundingVBMVsVoted_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVBMVsVoted_long$did <- dta_ReliefSurroundingVBMVsVoted_long$Relief * dta_ReliefSurroundingVBMVsVoted_long$d14

#table(dta_ReliefSurroundingVBMVsVoted_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVBMVsVotedSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVBMVsVoted_long)

# summary(m_ReliefVBMVsVotedSurrounding)

m_ReliefVBMVsVotedSurrounding <- coeftest(m_ReliefVBMVsVotedSurrounding, vcov=vcovHC(m_ReliefVBMVsVotedSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVBMVsVoted_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.83921 -0.29004 -0.06949  0.32171  1.13010 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.0642123  0.0014098  45.548  < 2e-16 ***
#   d14          0.4841377  0.0011497 421.089  < 2e-16 ***
#   Relief      -0.0308299  0.0012139 -25.398  < 2e-16 ***
#   did         -0.0894118  0.0017078 -52.354  < 2e-16 ***
#   AgeCat20142  0.0644657  0.0013081  49.282  < 2e-16 ***
#   AgeCat20143  0.1757983  0.0012392 141.860  < 2e-16 ***
#   AgeCat20144  0.2716895  0.0013413 202.563  < 2e-16 ***
#   gender2     -0.0175035  0.0008595 -20.365  < 2e-16 ***
#   Dem         -0.0283590  0.0010499 -27.011  < 2e-16 ***
#   NPA         -0.1230074  0.0012653 -97.216  < 2e-16 ***
#   Third       -0.0421348  0.0059824  -7.043 1.88e-12 ***
#   Black        0.0191739  0.0012221  15.689  < 2e-16 ***
#   Hispanic    -0.0190555  0.0027083  -7.036 1.98e-12 ***
#   Other       -0.0229688  0.0020014 -11.476  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4122 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.2718,	Adjusted R-squared:  0.2718 
# F-statistic: 2.699e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted by mail in 2012 were less likely to vote (by any method) in 2014 than were those living in counties not hit by michael who voted by mail in 2012


# VBM Only 2012 vs. VBM Only 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVBM2012VsVBM2014_long <- reshape(dta_ReliefSurrounding,
                                                      varying = c("VBMOnly2012", "VBMOnly2014"),
                                                      v.names = "voted",
                                                      timevar = "year",
                                                      times = c("VBMOnly2012", "VBMOnly2014"),
                                                      new.row.names = 1:946286,    #multiply the number of obs by 2
                                                      direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVBM2012VsVBM2014_long$year2014 <- ifelse(dta_ReliefSurroundingVBM2012VsVBM2014_long$year=='VBMOnly2014', 2014, 2012)
dta_ReliefSurroundingVBM2012VsVBM2014_long$d14 <- ifelse(dta_ReliefSurroundingVBM2012VsVBM2014_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVBM2012VsVBM2014_long <- data.table(dta_ReliefSurroundingVBM2012VsVBM2014_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVBM2012VsVBM2014_long$did <- dta_ReliefSurroundingVBM2012VsVBM2014_long$Relief * dta_ReliefSurroundingVBM2012VsVBM2014_long$d14

table(dta_ReliefSurroundingVBM2012VsVBM2014_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVBM2012VsVBM2014Surrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVBM2012VsVBM2014_long)

# summary(m_ReliefVBM2012VsVBM2014Surrounding)

m_ReliefVBM2012VsVBM2014Surrounding <- coeftest(m_ReliefVBM2012VsVBM2014Surrounding, vcov=vcovHC(m_ReliefVBM2012VsVBM2014Surrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVBM2012VsVBM2014_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.29267 -0.14745 -0.10070 -0.06078  1.00752 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.1082075  0.0011434  94.638  < 2e-16 ***
#   d14         -0.0178466  0.0009325 -19.139  < 2e-16 ***
#   Relief      -0.0188105  0.0009845 -19.106  < 2e-16 ***
#   did         -0.0186397  0.0013851 -13.457  < 2e-16 ***
#   AgeCat20142  0.0162704  0.0010609  15.336  < 2e-16 ***
#   AgeCat20143  0.0580571  0.0010051  57.764  < 2e-16 ***
#   AgeCat20144  0.1839781  0.0010878 169.127  < 2e-16 ***
#   gender2     -0.0185291  0.0006971 -26.581  < 2e-16 ***
#   Dem         -0.0049673  0.0008515  -5.834 5.43e-09 ***
#   NPA         -0.0377343  0.0010262 -36.771  < 2e-16 ***
#   Third       -0.0064438  0.0048520  -1.328 0.184151    
# Black       -0.0034370  0.0009912  -3.468 0.000525 ***
#   Hispanic    -0.0041676  0.0021965  -1.897 0.057778 .  
# Other        0.0004833  0.0016232   0.298 0.765907    
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3343 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.04755,	Adjusted R-squared:  0.04753 
# F-statistic:  3610 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted by mail in 2012 were less likely to vote by mail in 2014 than were those living in counties not hit by michael who voted by mail in 2012


# VBM Only 2012 vs. EIP Only 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingVBMVsEIP_long <- reshape(dta_ReliefSurrounding,
                                              varying = c("VBMOnly2012", "EIPOnly2014"),
                                              v.names = "voted",
                                              timevar = "year",
                                              times = c("VBMOnly2012", "EIPOnly2014"),
                                              new.row.names = 1:946286,    #multiply the number of obs by 2
                                              direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingVBMVsEIP_long$year2014 <- ifelse(dta_ReliefSurroundingVBMVsEIP_long$year=='EIPOnly2014', 2014, 2012)
dta_ReliefSurroundingVBMVsEIP_long$d14 <- ifelse(dta_ReliefSurroundingVBMVsEIP_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVBMVsEIP_long <- data.table(dta_ReliefSurroundingVBMVsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVBMVsEIP_long$did <- dta_ReliefSurroundingVBMVsEIP_long$Relief * dta_ReliefSurroundingVBMVsEIP_long$d14

# table(dta_ReliefSurroundingVBMVsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVBMVsEIPSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingVBMVsEIP_long)

# summary(m_ReliefVBMVsEIPSurrounding)

m_ReliefVBMVsEIPSurrounding <- coeftest(m_ReliefVBMVsEIPSurrounding, vcov=vcovHC(m_ReliefVBMVsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingVBMVsEIP_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.41301 -0.23972 -0.17069 -0.05028  1.01387 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.0986793  0.0013427  73.493  < 2e-16 ***
#   d14          0.0896726  0.0010950  81.891  < 2e-16 ***
#   Relief      -0.0220523  0.0011561 -19.074  < 2e-16 ***
#   did          0.0502690  0.0016266  30.905  < 2e-16 ***
#   AgeCat20142  0.0298631  0.0012459  23.970  < 2e-16 ***
#   AgeCat20143  0.1091343  0.0011803  92.465  < 2e-16 ***
#   AgeCat20144  0.1784731  0.0012774 139.711  < 2e-16 ***
#   gender2     -0.0096001  0.0008186 -11.727  < 2e-16 ***
#   Dem         -0.0260342  0.0009999 -26.036  < 2e-16 ***
#   NPA         -0.0741845  0.0012051 -61.559  < 2e-16 ***
#   Third       -0.0314731  0.0056978  -5.524 3.32e-08 ***
#   Black        0.0179732  0.0011640  15.441  < 2e-16 ***
#   Hispanic    -0.0045110  0.0025794  -1.749 0.080319 .  
# Other       -0.0067081  0.0019062  -3.519 0.000433 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3926 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.0549,	Adjusted R-squared:  0.05489 
# F-statistic:  4201 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted by mail in 2012 were MORE likely to vote early-in-person in 2014 than were those living in counties not hit by michael who voted by mail in 2012


# EIP Only 2012 vs. EIP Only 2014

# Reshape dataset from wide to long
dta_ReliefSurroundingEIPVsEIP_long <- reshape(dta_ReliefSurrounding,
                                              varying = c("EIPOnly2012", "EIPOnly2014"),
                                              v.names = "voted",
                                              timevar = "year",
                                              times = c("EIPOnly2012", "EIPOnly2014"),
                                              new.row.names = 1:946286,    #multiply the number of obs by 2
                                              direction = "long")



# Create year variable: 2014 = 1, 2012 = 0 
dta_ReliefSurroundingEIPVsEIP_long$year2014 <- ifelse(dta_ReliefSurroundingEIPVsEIP_long$year=='EIPOnly2014', 2014, 2012)
dta_ReliefSurroundingEIPVsEIP_long$d14 <- ifelse(dta_ReliefSurroundingEIPVsEIP_long$year2014==2014, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingEIPVsEIP_long <- data.table(dta_ReliefSurroundingEIPVsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingEIPVsEIP_long$did <- dta_ReliefSurroundingEIPVsEIP_long$Relief * dta_ReliefSurroundingEIPVsEIP_long$d14

# table(dta_ReliefSurroundingEIPVsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefEIPVsEIPSurrounding <- lm(voted ~ d14 + Relief + did + AgeCat2014 + gender2 + Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty + Holmes + Wakulla + Leon + Walton, data = dta_ReliefSurroundingEIPVsEIP_long)

# summary(m_ReliefEIPVsEIPSurrounding)

m_ReliefEIPVsEIPSurrounding <- coeftest(m_ReliefEIPVsEIPSurrounding, vcov=vcovHC(m_ReliefEIPVsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d14 + Relief + did + AgeCat2014 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other + Washington + Jackson + Calhoun + Gulf + Franklin + Gadsden + Liberty, data = dta_ReliefSurroundingEIPVsEIP_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.4379 -0.3322 -0.2242  0.5918  0.9436 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.2620006  0.0015282 171.441  < 2e-16 ***
#   d14         -0.0771095  0.0012463 -61.869  < 2e-16 ***
#   Relief      -0.0010990  0.0013159  -0.835    0.404    
# did          0.0274050  0.0018513  14.803  < 2e-16 ***
#   AgeCat20142  0.0392869  0.0014180  27.706  < 2e-16 ***
#   AgeCat20143  0.1522662  0.0013434 113.347  < 2e-16 ***
#   AgeCat20144  0.1555562  0.0014540 106.988  < 2e-16 ***
#   gender2     -0.0049238  0.0009317  -5.285 1.26e-07 ***
#   Dem         -0.0340373  0.0011381 -29.907  < 2e-16 ***
#   NPA         -0.1044403  0.0013716 -76.144  < 2e-16 ***
#   Third       -0.0505844  0.0064851  -7.800 6.19e-15 ***
#   Black        0.0202902  0.0013248  15.315  < 2e-16 ***
#   Hispanic     0.0032636  0.0029358   1.112    0.266    
# Other       -0.0190862  0.0021696  -8.797  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4468 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.03768,	Adjusted R-squared:  0.03767 
# F-statistic:  2832 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted early-in-person in 2012 were MORE likely to vote early-in-person in 2014 than were those living in counties not hit by michael who voted early-in-person in 2012

























# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)   # 90% CI
interval2 <- -qnorm((1-0.95)/2)  # 95% CI


# Put model estimates into temporary data.frames:
#================#
# Relief Policy  #
#================#

m_ReliefSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                         Coefficient = m_ReliefSurrounding[4, 1],
                                         SE = m_ReliefSurrounding[4, 2],
                                         modelName = "Voted 2012 vs. Voted 2014")

m_ReliefVBMVsVotedSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                   Coefficient = m_ReliefVBMVsVotedSurrounding[4, 1],
                                                   SE = m_ReliefVBMVsVotedSurrounding[4, 2],
                                                   modelName = "VBM 2012 vs. Voted 2014")

m_ReliefEDOnly2012VsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                        Coefficient = m_ReliefEDOnly2012VsEIPSurrounding[4, 1],
                                                        SE = m_ReliefEDOnly2012VsEIPSurrounding[4, 2],
                                                        modelName = "ED 2012 vs. EIP 2014")

# m_ReliefAnyEarlyVsVotedSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
#                                                       Coefficient = m_ReliefAnyEarlyVsVotedSurrounding[4, 1],
#                                                       SE = m_ReliefAnyEarlyVsVotedSurrounding[4, 2],
#                                                       modelName = "Any Early 2012 vs. Voted 2014")

m_ReliefEIPVsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                 Coefficient = m_ReliefEIPVsEIPSurrounding[4, 1],
                                                 SE = m_ReliefEIPVsEIPSurrounding[4, 2],
                                                 modelName = "EIP 2012 vs. EIP 2014")

m_ReliefVBM2012VsVBM2014Surrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                         Coefficient = m_ReliefVBM2012VsVBM2014Surrounding[4, 1],
                                                         SE = m_ReliefVBM2012VsVBM2014Surrounding[4, 2],
                                                         modelName = "VBM 2012 vs. VBM 2014")

m_ReliefVBMVsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                 Coefficient = m_ReliefVBMVsEIPSurrounding[4, 1],
                                                 SE = m_ReliefVBMVsEIPSurrounding[4, 2],
                                                 modelName = "VBM 2012 vs. EIP 2014")

m_ReliefVoted12VsAnyEarlySurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                          Coefficient = m_ReliefVoted12VsAnyEarlySurrounding[4, 1],
                                                          SE = m_ReliefVoted12VsAnyEarlySurrounding[4, 2],
                                                          modelName = "Voted 2012 vs. Any Early 2014")

m_ReliefVoted12VsEDSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                    Coefficient = m_ReliefVoted12VsEDSurrounding[4, 1],
                                                    SE = m_ReliefVoted12VsEDSurrounding[4, 2],
                                                    modelName = "Voted 2012 vs. ED 2014")

m_ReliefVoted12VsVBMSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                     Coefficient = m_ReliefVoted12VsVBMSurrounding[4, 1],
                                                     SE = m_ReliefVoted12VsVBMSurrounding[4, 2],
                                                     modelName = "Voted 2012 vs. VBM 2014")

m_ReliefVoted12VsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                     Coefficient = m_ReliefVoted2012VsEIPSurrounding[4, 1],
                                                     SE = m_ReliefVoted2012VsEIPSurrounding[4, 2],
                                                     modelName = "Voted 2012 vs. EIP 2014")

allModelFrame_ReliefSurrounding <- data.frame(rbind(m_ReliefEDOnly2012VsEIPSurrounding_frame1, m_ReliefEIPVsEIPSurrounding_frame1,
                                                    m_ReliefVBMVsEIPSurrounding_frame1, m_ReliefVBM2012VsVBM2014Surrounding_frame1,
                                                    m_ReliefVBMVsVotedSurrounding_frame1, m_ReliefVoted12VsAnyEarlySurrounding_frame1,
                                                    m_ReliefVoted12VsEDSurrounding_frame1, m_ReliefVoted12VsEIPSurrounding_frame1, 
                                                    m_ReliefVoted12VsVBMSurrounding_frame1, m_ReliefSurrounding_frame1))

# Plot
p_ReliefSurrounding <- ggplot(allModelFrame_ReliefSurrounding, aes(colour = Variable))
p_ReliefSurrounding <- p_ReliefSurrounding + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p_ReliefSurrounding <- p_ReliefSurrounding + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                                                ymax = Coefficient + SE*interval1),
                                                            lwd = 1, position = position_dodge(width = 1))

p_ReliefSurrounding <- p_ReliefSurrounding + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                                                 ymax = Coefficient + SE*interval2),
                                                             lwd = 0.5, position = position_dodge(width = 1),
                                                             shape = allModelFrame_ReliefSurrounding$Variable)
#p_ReliefSurrounding <- p_ReliefSurrounding + coord_flip() + theme_bw()
# p_ReliefSurrounding <- p_ReliefSurrounding + facet_grid(modelName ~. , rows = 4, cols = 2)
p_ReliefSurrounding <- p_ReliefSurrounding + facet_wrap(~modelName, nrow = 4)
p_ReliefSurrounding <- p_ReliefSurrounding + theme_bw(base_line_size = 0.2) + ggtitle("") + theme(plot.title = element_text(hjust = 0.5))
p_ReliefSurrounding <- p_ReliefSurrounding + theme(legend.title=element_blank(), 
                                                   axis.title.x=element_blank(),
                                                   axis.text.x=element_blank(),
                                                   axis.ticks.x=element_blank()) 
p_ReliefSurrounding <- p_ReliefSurrounding + scale_y_continuous() + ylab("")
p_ReliefSurrounding <- p_ReliefSurrounding + scale_color_manual(values=c('#000000','#000000', '#000000', '#000000')) 
p_ReliefSurrounding <- p_ReliefSurrounding + guides(colour = guide_legend(override.aes = list(shap_ReliefSurroundinge = c(1,2,3,4)),
                                                                          ncol = 2, nrow = 2))
p_ReliefSurrounding <- p_ReliefSurrounding + theme(legend.position='bottom', legend.box.background = element_rect(colour = "black"))
print(p_ReliefSurrounding)

ggsave("ReliefCountyFixedEffects.png", height=10, width=15, units='in', dpi=1000)






#=====================#
# Model Output/ LaTeX #  
#=====================#

library(lmtest) ; library(sandwich)
library(texreg)
library(xtable)


# Full DiD Models for Relief
ReliefModelsPart1 <- 
  texreg(list(m_ReliefSurrounding, m_ReliefEDOnly2012VsEIPSurrounding, m_ReliefEIPVsEIPSurrounding, m_ReliefVBM2012VsVBM2014Surrounding, m_ReliefVBMVsEIPSurrounding),
         caption="Relief Full DiD 2012-2014 Part 1",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("Voted 2012 Vs Voted 2014", "ED 2012 Vs EIP 2014", "EIP 2012 Vs EIP 2014", "VBM 2012 Vs VBM 2014", "VBM 2012 Vs EIP 2014"),
         override.se=list(summary(m_ReliefSurrounding)$coef[,2],
                          summary(m_ReliefEDOnly2012VsEIPSurrounding)$coef[,2],
                          summary(m_ReliefEIPVsEIPSurrounding)$coef[,2],
                          summary(m_ReliefVBM2012VsVBM2014Surrounding)$coef[,2],
                          summary(m_ReliefVBMVsEIPSurrounding)$coef[,2],
                          override.pval=list(summary(m_ReliefVBMVsEIPSurrounding)$coef[,4],
                                             summary(m_ReliefSurrounding)$coef[,4],
                                             summary(m_ReliefEDOnly2012VsEIPSurrounding)$coef[,4],
                                             summary(m_ReliefEIPVsEIPSurrounding)$coef[,4],
                                             summary(m_ReliefVBM2012VsVBM2014Surrounding)$coef[,4])))

ReliefModelsPart2 <- 
  texreg(list(m_ReliefVBMVsVotedSurrounding, m_ReliefVoted12VsAnyEarlySurrounding, m_ReliefVoted12VsEDSurrounding, m_ReliefVoted12VsVBMSurrounding, m_ReliefVoted2012VsEIPSurrounding),
         caption="Relief Full DiD 2012-2014 Part 2",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("VBM 2012 Vs Voted 2014", "Voted 2012 Vs Any Early 2014", "Voted 2012 Vs ED 2014", "Voted 2012 Vs VBM 2014", "Voted 2012 Vs EIP 2014"),
         override.se=list(summary(m_ReliefVBMVsVotedSurrounding)$coef[,2],
                          summary(m_ReliefVoted12VsAnyEarlySurrounding)$coef[,2],
                          summary(m_ReliefVoted12VsEDSurrounding)$coef[,2],
                          summary(m_ReliefVoted12VsVBMSurrounding)$coef[,2],
                          summary(m_ReliefVoted2012VsEIPSurrounding)$coef[,2],
                          override.pval=list(summary(m_ReliefVBMVsVotedSurrounding)$coef[,4],
                                             summary(m_ReliefVoted12VsAnyEarlySurrounding)$coef[,4],
                                             summary(m_ReliefVoted12VsEDSurrounding)$coef[,4],
                                             summary(m_ReliefVoted12VsVBMSurrounding)$coef[,4],
                                             summary(m_ReliefVoted2012VsEIPSurrounding)$coef[,4])))

# Full DiD Models for Relief (Robust Standard Errors)

ReliefModelsPart1 <- 
  texreg(list(m_ReliefSurrounding, m_ReliefEDOnly2012VsEIPSurrounding, m_ReliefEIPVsEIPSurrounding, m_ReliefVBM2012VsVBM2014Surrounding, m_ReliefVBMVsEIPSurrounding),
         caption="Relief Full DiD 2012-2014 Part 1",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("Voted 2012 Vs Voted 2014", "ED 2012 Vs EIP 2014", "EIP 2012 Vs EIP 2014", "VBM 2012 Vs VBM 2014", "VBM 2012 Vs EIP 2014"),
         override.se=list(m_ReliefSurrounding[,2],
                          m_ReliefEDOnly2012VsEIPSurrounding[,2],
                          m_ReliefEIPVsEIPSurrounding[,2],
                          m_ReliefVBM2012VsVBM2014Surrounding[,2],
                          m_ReliefVBMVsEIPSurrounding[,2],
                          override.pval=list(m_ReliefVBMVsEIPSurrounding[,4],
                                             m_ReliefSurrounding[,4],
                                             m_ReliefEDOnly2012VsEIPSurrounding[,4],
                                             m_ReliefEIPVsEIPSurrounding[,4],
                                             m_ReliefVBM2012VsVBM2014Surrounding[,4])))

ReliefModelsPart2 <- 
  texreg(list(m_ReliefVBMVsVotedSurrounding, m_ReliefVoted12VsAnyEarlySurrounding, m_ReliefVoted12VsEDSurrounding, m_ReliefVoted12VsVBMSurrounding, m_ReliefVoted2012VsEIPSurrounding),
         caption="Relief Full DiD 2012-2014 Part 2",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("VBM 2012 Vs Voted 2014", "Voted 2012 Vs Any Early 2014", "Voted 2012 Vs ED 2014", "Voted 2012 Vs VBM 2014", "Voted 2012 Vs EIP 2014"),
         override.se=list(m_ReliefVBMVsVotedSurrounding[,2],
                          m_ReliefVoted12VsAnyEarlySurrounding[,2],
                          m_ReliefVoted12VsEDSurrounding[,2],
                          m_ReliefVoted12VsVBMSurrounding[,2],
                          m_ReliefVoted2012VsEIPSurrounding[,2],
                          override.pval=list(m_ReliefVBMVsVotedSurrounding[,4],
                                             m_ReliefVoted12VsAnyEarlySurrounding[,4],
                                             m_ReliefVoted12VsEDSurrounding[,4],
                                             m_ReliefVoted12VsVBMSurrounding[,4],
                                             m_ReliefVoted2012VsEIPSurrounding[,4])))

